\chapter{Analisi dei risultati}

In questo capitolo mostriamo i risultati numerici ottenuti risolvendo problemi test tramite l'utilizzo della libreria \texttt{Heat-Spectral}.
Esiste un'enorme varietà di problemi che possono essere simulati per mezzo del nostro codice, a patto che il dominio consista in una regione rettangolare.
Cominciamo, pertanto, con l'illustrare i risultati relativi a un esempio piuttosto semplice, di cui conosciamo la soluzione analitica, in modo da mostrare le proprietà di convergenza dei metodi numerici utilizzati.
In un secondo momento, poi, passeremo a descrivere la simulazione di un problema concreto, relativo all'ambito ingegneristico della trasmissione di calore.
Naturalmente, in questo caso, non abbiamo a disposizione la soluzione analitica e, pertanto, intendiamo unicamente mostrare le potenzialità della libreria che abbiamo costruito, illustrando in che modo sia possibile implementare problemi più sofisticati.


\section{Caso test con soluzione esatta}

Come primo passo è stato necessario costruire un problema parabolico corredato dalla sua soluzione esatta, in modo da poter sviscerare alcune considerazioni relative ai risultati di convergenza.
Dunque, cominciamo con lo scegliere un dominio quadrato $\varOmega=[0,2\pi]\times[0,2\pi]$ sul quale, per semplicità, cerchiamo di risolvere per via analitica il problema differenziale
\[
	u_t - \varDelta u = 0,
\]
corredato dalla condizione iniziale e da condizioni al bordo di tipo Dirichlet omogenee.
Scegliamo di risolvere il problema analitico attraverso il metodo di separazione di variabili.
In questo modo sappiamo che la soluzione potrà essere espressa in serie, secondo la formula
\[
	u(\vec{x},t) = \sum_{k=0}^\infty c_k(t) u_k(\vec{x}).
\]
Ora, considerando quest'ultima espressione come soluzione del problema, inseriamola nell'equazione differenziale in modo che, per ogni $k$, otteniamo
\[
	c'_k(t)u_k(\vec{x})-c_k(t)\varDelta u_k(\vec{x}) = 0.
\]
A questo punto, raggruppando i termini relativi alla medesima variabile, ricaviamo la relazione
\[
	-\dfrac{c'_k(t)}{c_k(t)} = -\dfrac{\varDelta u_k(\vec{x})}{u_k(\vec{x})} = \lambda_k
\]
e i primi due termini devono essere uguali a una costante proprio perché coinvolgono variabili diverse.
La prima equazione è di risoluzione immediata: si ottiene che
\[
	c_k(t) = e^{-\lambda_k t},
\]
a meno di una costante moltiplicativa, da fissare per mezzo della condizione iniziale.
Invece, per quanto riguarda la seconda, otteniamo
\[
	-\varDelta u_k(\vec{x}) = \lambda_k u_k(\vec{x}),
\]
cioè l'equazione per gli autovalori e gli autovettori del laplaciano, con condizioni al bordo di tipo Dirichlet omogenee.
Grazie al metodo di separazione di variabili, e all'imposizione delle condizioni al contorno, ricaviamo immediatamente che le autofunzioni devono essere del tipo
\[
	u_{n,m}(\vec{x}) = \sin(nx)\sin(my),\qquad n,m\in\N
\]
e, di conseguenza, l'autovalore associato è definito come
\[
	\lambda_{n,m} = n^2+m^2,\qquad n,m\in\N.
\]
Dopo tutte queste considerazioni, siamo quasi pronti a fabbricare il nostro caso d'esempio.
Per semplificare ulteriormente il problema, scegliamo come dato iniziale proprio un'autofunzione.
Quindi consideriamo
\[
	u(\vec{x},0) = \sin(nx)\sin(my).
\]
In questo modo, grazie ai conti svolti in precedenza, riusciamo a ottenere la soluzione esatta del problema differenziale, che risulta essere
\[
	u(\vec{x},t) = e^{-(n^2+m^2)t}\sin(nx)\sin(my).
\]
A questo punto siamo finalmente pronti per utilizzare il nostro codice e risolvere il problema dal punto di vista numerico, in modo da avere tutte le informazioni necessarie per le analisi successive.\\
Assegniamo dei valori ai parametri $n$ ed $m$, in modo da selezionare un problema ben preciso.
Optiamo per la scelta $n=m=2$, in modo da ottenere un problema evolutivo che non converga troppo velocemente verso la soluzione asintotica nulla.
Di fatto, il problema non è nient'altro che l'equazione del calore, già analizzata nel capitolo precedente:
\begin{equation}
	\left\{
	\begin{array}{ll}
		u_t - \varDelta u = 0		&\qquad\textrm{in}\ \varOmega \times [0,T]\\[2mm]
		u(x,y,0) = \sin(2x) \sin(2y)	&\qquad\textrm{in}\ \varOmega\\[2mm]
		u(x,y,t) = 0				&\qquad\textrm{su}\ \partial\varOmega,
	\end{array}
	\right.
\end{equation}
la cui soluzione esatta è, come ricavato in precedenza,
\[
	u(x,y,t) = e^{-8t} \sin(2x) \sin(2y).
\]

\subsection{Definizione del problema variazionale}

Dal punto di vista del codice di calcolo, la definizione del problema è collocata nel file \texttt{exact.hpp}.
Innanzitutto vengono definite la soluzione esatta, le sue derivate spaziali e il funzionale, in modo da poter calcolare le norme degli errori.
\begin{c++}
struct ExactSolution {
    double operator()(double x,double y,double t) const {
        return std::exp(-8*t)*std::sin(2*x)*std::sin(2*y);
    }
};

struct ExactSolutionX {
    double operator()(double x,double y,double t) const {
        return 2*std::exp(-8*t)*std::cos(2*x)*std::sin(2*y);
    }
};

struct ExactSolutionY {
    double operator()(double x,double y,double t) const {
        return 2*std::exp(-8*t)*std::sin(2*x)*std::cos(2*y);
    }
};

struct ExactFunc {
    double operator()(double x,double y) const {
        return 0;
    }
};
\end{c++}
Una volta definiti questi elementi, viene costruita la forma bilineare, tenendo anche in considerazione la condizione iniziale e le condizioni al contorno:
\begin{c++}
template<class MESH>
class ExactBilinear {
    public:
        ExactBilinear(MESH const & me): me(me) {
            for(int i=0;i<4;++i) {
                bcond[i]=spectral::dirichlet;
            }
        }

        template<class U,class V> inline
        double operator()(U const & u,V const & v,int n) const {
            return
                + me.prodL2(me.Dx(u,n),me.Dx(v,n),n)
                + me.prodL2(me.Dy(u,n),me.Dy(v,n),n);
        }

        // Initial data
        struct initial{
            double operator()(double const & x,
                              double const & y) const {
                return std::sin(2*x)*std::sin(2*y);
            }
        };
        
        // Boundary conditions
        spectral::bctype bcond[4];

        // Function for the Dirichlet boundary condition
        struct dirichlet{
            double operator()(double const & x,
	                      double const & y) const {
                return 0;
            }
        };

        // Method on
        template<unsigned int N, typename T=void>
        struct on : public dirichlet {};

    private:
        MESH const & me;
};
\end{c++}
È bene sottolineare che, in un problema elementare come quello preso in considerazione, la struttura preposta a gestire le condizioni al contorno risulta fin troppo elaborata.
Di certo, nel caso di problemi semplici, sarebbe possibile trovare un modo più intuitivo di gestire queste condizioni, senza dover necessariamente ricorrere all'impiego dei template.
In realtà, però, il nostro codice ha la pretesa di gestire contemporaneamente condizioni essenziali e naturali, eventualmente imponendone di un certo tipo su alcuni bordi e di tipo diverso su altri.
In questo modo, per avere la possibilità di usufruire di potenzialità più sofisticate, è naturale dover complicare lo stile di programmazione, in modo da ottenere un codice che risulti organicamente coerente.
Nel caso in cui si abbia la necessità di imporre sia condizioni di Dirichlet sia di Neumann, bisogna selezionare il tipo corretto all'interno del vettore \texttt{bcond}.
Dopodiché si creano diversi metodi \texttt{on}, che richiamano le funzioni corrispondenti al dato essenziale o naturale, templetizzate sul lato del dominio a cui fanno riferimento.\\
Proseguendo la descrizione del nostro esempio, infine, mostriamo come il dominio spaziale sia descritto semplicemente definendone due estremi: il vertice in basso a destra e quello in alto a sinistra.
\begin{c++}
const spectral::Point<double> LA(0.0,0.0);
const spectral::Point<double> LB(2.0*M_PI,2.0*M_PI);
\end{c++}

\subsection{Scelta dei parametri e risoluzione numerica}

Una volta compiute tutte queste operazioni preliminari, finalmente giungiamo alla costruzione del \texttt{main()}, che sarà contenuto nel file \texttt{exact.cpp}.
Per prima cosa è necessario assegnare i valori numerici ai parametri della classe \texttt{Problem}.
\begin{c++}
// Number of elements along the x-axis
int nx(2);
// Number of elements along the y-axis
int ny(2);
// Number of nodes
int nnodes(4);
// Time step
double deltat(0.05);
// Number of time steps
int ntime(20);
// Relaxation parameter
double theta(0.5);
// nmax for the newton method
int newton_nmax((int)1e4);
// Order of the problem in the week form
int order(1);
// Tollerance paramethers for the newton method
double legDelta(1e-15);
double legEps(1e-15);
\end{c++}
In questo particolare caso abbiamo assegnato $2$ elementi spettrali su ciascun asse, $4$ nodi di interpolazione per i polinomi su ciascun elemento, un valore per $\Delta t$ pari a $0.05$ secondi e $20$ passi temporali, in modo che l'istante finale della simulazione corrisponda a $1$ secondo.
Inoltre, abbiamo imposto $\vartheta=0.5$, in modo da utilizzare il metodo di Crank-Nicolson.
Infine, i restanti parametri servono a regolare le tolleranze nel calcolo delle formule di quadratura.\\
A questo punto è possibile procedere con la risoluzione vera e propria del problema, ma non prima di aver definito le variabili locali necessarie alla memorizzazione dei risultati.
Tralasciamo questi dettagli e passiamo a mostrare come si inizializza il problema, oltre a descrivere in che modo sia possibile calcolare gli errori nelle diverse norme.
\begin{c++}
Problem<EXACT_PROBLEM> pb(LA, LB, nx, ny, nnodes, deltat, theta,
                          newton_nmax, legDelta, legEps, order);
pb.initialize();
L2error = sqrt(max(0.,pb.evaluate<myL2error>(0)));
H1error = sqrt(max(0.,pb.evaluate<myH1error>(0)));
LinfL2error = L2error;
L2H1error = deltat/2 * pow(H1error,2);
l2errorfile << L2error << endl;
h1errorfile << H1error << endl;
\end{c++}
Gli errori vengono calcolati dal metodo \texttt{evaluate}, che prende in ingresso l'istante temporale a cui la risoluzione del problema è giunta, oltre a un parametro come template, che si riferisce al tipo di errore che si vuole ottenere.
Il metodo restituisce i valori corretti ma elevati al quadrato dunque, per ottenere il risultato cercato, è necessario estrarre una radice, con l'accorgimento di annullare gli eventuali valori negativi che potrebbero risultare da errori numerici.
Notiamo che terremo traccia degli errori ad ogni passo temporale, in modo da poter poi svolgere alcune analisi sugli andamenti di questi valori.
A tal proposito è necessario utilizzare gli \texttt{ofstream} per memorizzare in alcuni file l'entità degli errori sopra citati.\\
Oltre a memorizzare gli errori di discretizzazione dei singoli problemi ellittici, che vengono risolti a ogni istante temporale, abbiamo intenzione di ricavare anche gli errori nelle norme \emph{classiche} utilizzate nei problemi parabolici: la norma $L^\infty([0,T],L^2(\varOmega))$ e la norma $L^2([0,T],H^1(\varOmega))$.
Per il calcolo di questi due valori è necessario, nel primo caso, tenere traccia del valore più grande su tutto l'intervallo temporale mentre, nel secondo, bisogna calcolare un integrale numericamente, attraverso un'opportuna formula di quadratura.
Dato che in questo caso la griglia di discretizzazione temporale è imposta dal metodo delle differenze finite, non abbiamo la possibilità di utilizzare formule di quadratura di grado elevato, poiché avremmo bisogno di molti nodi in più rispetto a quelli che abbiamo realmente a disposizione.
In questo caso, il metodo dei trapezi sembra essere un buon compromesso, poiché si tratta di una formula di quadratura di grado sufficientemente accurato e, contemporaneamente, di facile implementazione.
Infatti, avendo a disposizione l'errore $e(t) = \Vert u(t) - u^*_N(t) \Vert_{H^1(\varOmega)}$ a ogni istante temporale, utilizziamo la formula
\[
	\int_0^T e(t)\,\D t \approx \dfrac{\Delta t}{2}(e_0 + e_T) + \Delta t \sum_{k=1}^{T-1} e_k,
\]
dove $\vec{e} = (e_0,\ldots,e_T) \in \R^{T+1}$ è la rappresentazione dell'errore $e(t)$ sui nodi di discretizzazione del $\vartheta$-metodo.
I termini dell'integrale numerico sono facilmente calcolabili all'interno del ciclo di avanzamento in tempo, facendo particolare attenzione ai valori estremi, per i quali il coefficiente è diverso.
Per questo, come mostrato nel segmento di codice, l'inizializzazione dell'errore in norma $L^2([0,T],H^1(\varOmega))$ utilizza una formula diversa da quella che figurerà all'interno del ciclo per la risoluzione in tempo.\\
Una volta svolte queste operazioni preliminari, possiamo finalmente risolvere il problema evolutivo, tralasciando nella descrizione tutte i metodi che consentono di salvare in vari file la soluzione, oltre a quelle che ne permettono la rappresentazione grafica.
\begin{c++}
for (int t=0; t<ntime; ++t) {
    actualtime = deltat*(t+1);
    pb.solve();
    // Computes the error in L2 norm:
    L2error = sqrt(max(0.,pb.evaluate<myL2error>(actualtime)));
    LinfL2error = max(LinfL2error,L2error);
    l2errorfile << L2error << endl;
    // Computes the error in H1 norm:
    H1error = sqrt(max(0.,pb.evaluate<myH1error>(actualtime)));
    L2H1error = L2H1error + (deltat * pow(H1error,2));
    h1errorfile << H1error << endl;
}
\end{c++}
Per evitare di inserire una condizione di tipo \texttt{if} all'interno del ciclo \texttt{for}, in modo da non gravare sull'onerosità computazionale, fino a questo punto non si è fatto caso al fatto che l'ultimo termine nella formula di quadratura per l'errore $L^2([0,T],H^1(\varOmega))$ ha un coefficiente dimezzato.
Il problema è risolto fuori dal ciclo, grazie all'istruzione
\begin{c++}
	L2H1error = L2H1error - (deltat/2 * pow(H1error,2));
\end{c++}
che rimette le cose a posto.
Di fatto, in questo modo, semplicemente eliminiamo il termine in eccesso, grazie al fatto che la variabile \texttt{H1error} non viene più modificata una volta usciti dal ciclo \texttt{for} in tempo.\\
Fin a qui ci siamo occupati di mostrare in che modo sia possibile impostare e risolvere numericamente un problema differenziale sfruttando la nostra libreria.
Oltre a questo, però, grazie ai metodi di \emph{post-processing} è immediato produrre dei grafici della soluzione, sia di tipo tridimensionale sia bidimensionali, ottenuti mappando la funzione in un piano.
Inoltre, abbiamo anche la possibilità di mostrare il \emph{pattern} della matrice, relativa al problema algebrico, al variare del numero di elementi spettrali.
In Figura~\ref{fig:pattern} possiamo confrontare la quantità di elementi non nulli che si ottengono, man mano che si infittisce la griglia spaziale, per esempio aumentando il numero di elementi.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Griglia di $2\times 2$ elementi}.]
		{\includegraphics[width=.45\columnwidth]{pattern2x2}} \quad
	\subfigure[\emph{Griglia di $10\times 10$ elementi}.]
		{\includegraphics[width=.45\columnwidth]{pattern10x10}}
	\caption{Immagini relative al \emph{pattern} delle matrici, per il problema algebrico derivato dal metodo SEM-NI.}
	\label{fig:pattern}
\end{figure}
Immagini di questo tipo si possono ottenere facilmente attraverso il metodo \texttt{printPattern} della nostra libreria, una volta accoppiato, ad esempio, con la funzione \texttt{spy} di Matlab oppure di Octave.
Nel nostro caso abbiamo scelto $4$ nodi di quadratura per i polinomi che compongono la base spettrale e, nell'immagine a risoluzione più bassa, questo fatto si può constatare in maniera diretta.\\
Naturalmente, a parte queste rappresentazioni che sono comunque relative a un'analisi più marginale, dicevamo che i metodi di \emph{post-processing} contenuti nel nostro codice consentono anche di ottenere dei grafici della soluzione a vari istanti temporali. 
Ad esempio, per quanto riguarda il problema che stiamo analizzando, siamo in grado di produrre le immagini di Figura~\ref{fig:iniziale} relative al dato iniziale del problema.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Rappresentazione bidimensionale}.]
		{\includegraphics[width=.45\columnwidth]{init2d}} \quad
	\subfigure[\emph{Rappresentazione tridimensionale}.]
		{\includegraphics[width=.45\columnwidth]{init3d}}
	\caption{Immagini relative al dato iniziale nel problema di esempio.}
	\label{fig:iniziale}
\end{figure}
In questo particolare caso abbiamo modificato la griglia, rispetto a quella definita in precedenza, rendendola più fitta.
In questo in modo riusciamo a ottenere delle immagini decisamente ben definite mentre, se avessimo mantenuto una griglia lasca, saremmo ugualmente andati a convergenza, ma la qualità delle immagini ne avrebbe chiaramente risentito.\\
La simulazione è stata protratta per un lasso di tempo di un secondo, poiché il problema è fortemente dissipativo e la soluzione converge a $0$ molto rapidamente.
Quindi, essendo chiaro che la soluzione al tempo finale sarà identicamente nulla, scegliamo di riportare in Figura~\ref{fig:after} la superficie che si ottiene dopo aver risolto il problema evolutivo per un decimo di secondo.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Rappresentazione bidimensionale}.]
		{\includegraphics[width=.45\columnwidth]{after2d}} \quad
	\subfigure[\emph{Rappresentazione tridimensionale}.]
		{\includegraphics[width=.45\columnwidth]{after3d}}
	\caption{Immagine della soluzione che si ottiene dopo $0.1$~secondi.}
	\label{fig:after}
\end{figure}
Come possiamo notare dalle immagini, il problema test che stiamo analizzando è decisamente banale, poiché la soluzione consiste unicamente nello smorzamento del dato iniziale, in modo che converga a un risultato nullo.
Di fatto è stata scelta questa equazione poiché, non essendo troppo complicata, ci ha permesso di calcolarne la soluzione esatta per via analitica, in modo da poter svolgere gli studi di convergenza presenti nella prossima sezione.

\subsection{Analisi di convergenza degli schemi numerici}

Una volta risolto il problema dal punto di vista numerico, e grazie alla soluzione esatta calcolata per via analitica, possiamo ora verificare che i risultati e le prestazioni del codice siano in linea con quanto previsto dalla teoria.
A questo scopo analizziamo l'andamento nel tempo delle norme $H^1(\varOmega)$ e $L^2(\varOmega)$ dell'errore tra soluzione numerica e soluzione esatta.
Le simulazioni dalle quali abbiamo ricavato le informazioni necessarie a svolgere questi conti sono state effettuate su una griglia di $10$ elementi per lato.
Inoltre, abbiamo cominciato le analisi con il metodo di Crank-Nicolson, suddividendo il dominio temporale di un secondo in $1000$ passi.
Impostando i parametri nel suddetto modo, riusciamo a ricavare i valori delle norme dell'errore nei singoli istanti temporali e a rappresentarli nei grafici di Figura~\ref{fig:timeCN}.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Errore in norma $L^2(\varOmega)$}.]
		{\includegraphics[width=.45\columnwidth]{L2timeCN}} \quad
	\subfigure[\emph{Errore in norma $H^1(\varOmega)$}.]
		{\includegraphics[width=.45\columnwidth]{H1timeCN}}
	\caption{Andamento delle norme dell'errore fra soluzione numerica e soluzione analitica a ogni istante temporale, quando scegliamo come risolutore il metodo di Crank-Nicolson.}
	\label{fig:timeCN}
\end{figure}
Rispetto a queste immagini ci sono diverse considerazioni da fare.
Dovendo procedere con ordine, cominciamo a discutere gli andamenti delle curve continue tralasciando, almeno per il momento, la descrizione e il significato di quelle tratteggiate.\\
La norma $L^2(\varOmega)$ dell'errore è ovviamente nulla nell'istante iniziale, dato che la soluzione numerica viene ottenuta attraverso un'interpolazione dalla condizione iniziale imposta.
Poi, risolvendo il problema nei successivi istanti temporali, la norma cresce fino a un picco, per poi decrescere in maniera esponenziale verso zero, indicando così la convergenza dello schema numerico.
Per quanto riguarda la norma $H^1(\varOmega)$, invece, essa assume il valore massimo nell'istante iniziale, dato che non consiste in una semplice valutazione di un integrale, ma contiene in sé anche il termine relativo al gradiente.
L'andamento esponenziale decrescente, tuttavia, è invece analogo al caso precedente.\\
L'aspetto di questi grafici non varia qualitativamente se modifichiamo i parametri del problema, relativi al passo spaziale e temporale, oppure se utilizziamo un $\vartheta$-metodo diverso, ovviamente nel caso in cui si garantisca la convergenza della soluzione numerica a quella esatta.
Per esempio, possiamo provare a ripetere l'analisi precedente, impiegando però il metodo di Eulero Implicito.
In questo caso riusciamo a ottenere i grafici degli errori rappresentati in Figura~\ref{fig:timeEI}.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Errore in norma $L^2(\varOmega)$}.]
		{\includegraphics[width=.45\columnwidth]{L2timeEI}} \quad
	\subfigure[\emph{Errore in norma $H^1(\varOmega)$}.]
		{\includegraphics[width=.45\columnwidth]{H1timeEI}}
	\caption{Andamento delle norme dell'errore fra soluzione numerica e soluzione analitica a ogni istante temporale, quando scegliamo come risolutore il metodo di Eulero Implicito.}
	\label{fig:timeEI}
\end{figure}
Osservando le immagini, si può immediatamente constatare quando accennato in precedenza, e cioè che le curve degli errori hanno grosso modo gli stessi intervalli di monotonia, a prescindere dallo schema numerico utilizzato.\\
Passiamo ora a discutere le differenze fra le diverse curve nello stesso grafico, al variare del colore o del tipo di tratto.
Nei grafici di Figura~\ref{fig:timeCN} e di Figura~\ref{fig:timeEI} sono rappresentate, con colori diversi, le curve degli errori corrispondenti a diverse scelte del passo temporale.
Passando dal blu al rosso, e infine al verde, stiamo dimezzando il passo di avanzamento in tempo, in modo da verificare l'ordine di convergenza dei vari schemi.
In particolare, grazie a queste immagini, possiamo mostrare come l'ordine di convergenza, coerentemente con quanto predetto dalla teoria, sia quadratico nel caso di Crank-Nicolson e lineare quando si sceglie di applicare il metodo di Eulero Implicito.
A questo punto non ci resta altro che descrivere il significato delle curve con un tratto diverso, contenute nelle immagini che stiamo analizzando.
Nel caso di Crank-Nicolson, le linee tratteggiate rappresentano la divisione per quattro della corrispondente curva con linea piena.
In questo modo si riesce a mostrare che, escludendo un breve transitorio, si ha una precisa corrispondenza fra le curve con passi temporali differenti.
Di fatto, il dimezzamento del passo temporale si riflette nella riduzione a un quarto della norma dell'errore, indicando perciò la convergenza quadratica dello schema numerico.
Lo stesso discorso vale per il metodo di Eulero Implicito, dove però le linee tratteggiate si riferiscono al dimezzamento del valore corrispondente alla linea piena.
A differenza di quanto accadeva nel caso di Crank-Nicolson, che è un metodo del secondo ordine, ora al dimezzamento del passo temporale corrisponde un dimezzamento della norma dell'errore.
Dunque, anche in questo caso ritroviamo un aspetto previsto dalla teoria, cioè che il metodo di Eulero Implicito garantisce una convergenza di tipo lineare.

Fin qui abbiamo analizzato il comportamento delle norme, negli spazi $L^2(\varOmega)$ e $H^1(\varOmega)$, relative all'errore di discretizzazione dei singoli problemi ellittici che vengono risolti a ogni istante temporale.
In realtà, come accennavamo in precedenza, le ambientazioni funzionali corrette per la soluzione numerica di un problema parabolico sono gli spazi $L^\infty([0,T],L^2(\varOmega))$ e $L^2([0,T],H^1(\varOmega))$.
Per questo motivo, nel seguito dell'analisi, mostreremo i risultati di convergenza relativi a queste norme, le quali tengono traccia dell'errore di discretizzazione sia in spazio sia in tempo.
A questo proposito ricordiamo che, poiché nei problemi evolutivi l'errore di discretizzazione è composto da entrambi i contributi, nell'analizzarne il comportamento sarà necessario isolare il fenomeno oggetto di studio, cercando di ridurre gli errori dovuti alla componente di disturbo.
Per esempio, relativamente allo studio di convergenza in funzione del passo della griglia spaziale, cercheremo di ridurre ai minimi termini l'errore dovuto alla discretizzazione con differenze finite.
Negli altri casi ci comporteremo in maniera analoga.\\
Venendo ora all'analisi vera e propria, decidiamo che i parametri rispetto ai quali mostreremo il comportamente dell'errore sono:
\begin{itemize}
	\item il numero di elementi spettrali in cui suddividere il dominio del problema;
	\item il numero di nodi di quadratura utilizzati nel calcolo dei polinomi di Legendre;
	\item l'ampiezza del passo temporale negli schemi alle differenze finite.
\end{itemize}
Cominciamo, pertanto, ad affrontare il primo problema: la discussione dei risultati di convergenza dello schema numerico man mano che il passo di griglia spaziale viene dimezzato.
Come già accennato, visto che per ora intendiamo occuparci esclusivamente del contributo dell'errore dovuto alla discretizzazione spaziale, scegliamo un passo di avanzamento temporale decisamente piccolo, $\Delta t = 0.001$, in modo da ridurre la componente di errore causata dal $\vartheta$-metodo.
Inoltre, scegliamo anche un grado polinomiale non troppo elevato, in modo da non rendere lo schema troppo convergente, poiché in quel caso non saremmo in grado di leggere degli errori sufficientemente elevati da permettere l'analisi.
Per questo motivo, relativamente a questo punto dello studio, limitiamo il numero di nodi per la formula di quadratura dei polinomi a $4$.
Dopo aver prestato attenzione a tutti questi accorgimenti, possiamo ottenere i grafici in Figura~\ref{fig:space}, riportanti i valori delle norme dell'errore al variare della griglia di discretizzazione spaziale.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Metodo di Crank-Nicolson}.]
		{\includegraphics[width=.45\columnwidth]{spaceCN}} \quad
	\subfigure[\emph{Metodo di Eulero Implicito}.]
		{\includegraphics[width=.45\columnwidth]{spaceEI}}
	\caption{Logaritmi degli errori di discretizzazione in funzione del logaritmo del passo della griglia spaziale. In blu riportiamo i valori calcolati con la norma $L^\infty([0,T],L^2(\varOmega))$, mentre in rosso quelli relativi alla norma $L^2([0,T],H^1(\varOmega))$.}
	\label{fig:space}
\end{figure}
Abbiamo preso in considerazione i casi con 2, 4, 8 e 16 elementi spettrali per lato e, dato che ciascun lato del dominio è lungo $2\pi$, ciò significa che i passi di discretizzazione studiati saranno $\pi$, $\pi/2$, $\pi/4$ e $\pi/8$.
I valori numerici a cui corrispondono i punti dei grafici sono riassunti in Tabella~\ref{tab:passo_spaziale}.
\begin{table}[tb]
	\centering
	\begin{tabular}{*4{c}}
		\toprule
		Metodo					& Passo		& Errore $L^\infty([0,T],L^2(\varOmega))$	& Errore $L^2([0,T],H^1(\varOmega))$\\
		\midrule
		\multirow{4}*{Crank-Nicolson}	& $3.141593$	& $5.12828\cdot 10^{-1}$				& $1.91763$\\
								& $1.570796$	& $1.47338\cdot 10^{-2}$				& $3.98731\cdot 10^{-1}$\\
								& $0.785398$	& $1.22036\cdot 10^{-3}$				& $4.06481\cdot 10^{-2}$\\
								& $0.392699$	& $5.15912\cdot 10^{-5}$				& $5.15587\cdot 10^{-3}$\\
		\midrule
		\multirow{4}*{Eulero Implicito}	& $3.141593$	& $5.18460\cdot 10^{-1}$				& $1.92410$\\
								& $1.570796$	& $1.01953\cdot 10^{-2}$				& $3.98285\cdot 10^{-1}$\\
								& $0.785398$	& $4.73587\cdot 10^{-3}$				& $4.12320\cdot 10^{-2}$\\
								& $0.392699$	& $4.60851\cdot 10^{-3}$				& $8.41726\cdot 10^{-3}$\\
		\bottomrule
	\end{tabular}
	\caption{Norme degli errori al variare del passo di discretizzazione spaziale.}
	\label{tab:passo_spaziale}
\end{table}
In realtà, i valori rappresentati nei grafici sono stati filtrati da una trasformazione logaritmica, applicata sia ai passi di discretizzazione, sia ai valori dei corrispondenti errori.
In questo modo, avendo trasformato entrambi gli assi in scala logaritmica, le curve degli errori in funzione del passo di griglia devono necessariamente essere delle rette, in accordo con quanto predetto dalla teoria.
Anche solo attraverso un rapido sguardo alla reciproca posizione dei vari punti, ci si può immediatamente accorgere che i nostri risultati rispecchiano senza alcun dubbio l'inquadramento teorico.
Unica eccezione sembra essere il valore corrispondente al passo di griglia più piccolo, nel caso di errore $L^\infty([0,T],L^2(\varOmega))$ calcolato con il metodo di Eulero Implicito.
In quel caso sembra che la componente di errore dovuta ad altri fenomeni sia decisamente più influente, poiché il valore non è in accordo con l'andamento lineare.
In ogni caso, per determinare dal punto di vista numerico gli ordini di convergenza, abbiamo effettuato una regressione lineare con l'ausilio del software scientifico \texttt{R}, ottimo strumento per indagare questioni di tipo statistico.
In questo modo siamo riusciti a ricavare che la distribuzione dei valori è effettivamente disposta lungo una retta, poiché il valore di $R^2$ per la regressione è prossimo a $1$ e tutte le ipotesi necessarie sono soddisfatte, garantendo la significatività di tale operazione.
L'unico caso dubbio è quello citato poco sopra, per il quale abbiamo scelto di calcolare la retta senza includere il valore estraneo al \emph{trend} lineare.
Grazie a questo strumento statistico siamo in grado di calcolare i coefficienti angolare delle varie rette di regressione, in modo da stabilire gli ordini di convergenza nelle diverse norme.
Per quanto riguarda la norma $L^\infty([0,T],L^2(\varOmega))$, ricaviamo un valore pari a $4.34309$ nel caso di Crank-Nicolson e $3.38723$ con Eulero Implicito.
Relativamente agli errori in norma $L^2([0,T],H^1(\varOmega))$, invece, gli ordini di convergenza sono pressoché simili nei due metodi: $2.89108$ con Crank-Nicolson e $2.67818$ con Eulero Implicito.\\
Abbiamo scelto di non effettuare alcuna analisi sul metodo di Eulero Esplicito poiché, considerate le ulteriori condizioni necessarie a garantirne la convergenza, non avremmo potuto effettuare lo studio con gli stessi parametri impiegati negli altri casi.
Ad ogni modo non riteniamo che questo comporti una limitazione nel lavoro, poiché abbiamo focalizzato la nostra attenzioni sugli altri studi di convergenza, i quali portano alla luce dei fenomeni decisamente più interessanti.

Un'analisi simile a quanto fatto fin'ora può essere condotta variando il grado dei polinomi dell'interpolazione numerica, e osservando in che modo cambiano le norme degli errori in questo scenario.
In realtà noi modificheremo opportunamente il numero di nodi di quadratura che vengono utilizzati per il calcolo dei polinomi e, di conseguenza, in questo modo riusciamo a cambiare il grado dei polinomi stessi.
In questo caso fissiamo a $2$ il numero di elementi spettrali, per una ragione analoga alla precedente.
Non vogliamo che una griglia fitta porti a $0$ gli errori di discretizzazione troppo in fretta, mascherando così il contributo alla riduzione dell'errore dovuto all'aumento del grado polinomiale.
Inoltre, proprio come nel caso precedente, manteniamo un passo di avanzamento temporale molto piccolo, sempre $\Delta t = 0.001$, in modo da minimizzare il contributo di errore in tempo.
Per questo tipo di analisi non riteniamo che il tipo di $\vartheta$-metodo utilizzato possa influire in modo significativo sui risultati.
Pertanto, limitatamente a questo studio, scegliamo di analizzare esclusivamente i dati ottenuti attraverso il metodo di Crank-Nicolson.
In queso modo, effettuando diverse simulazioni per ottenere tutti le informazioni necessarie all'analisi, riusciamo a ricavare i valori che abbiamo inserito in Tabella~\ref{tab:grado}.
\begin{table}[tb]
	\centering
	\begin{tabular}{*3{c}}
		\toprule
		Nodi per i polinomi	& Errore $L^\infty([0,T],L^2(\varOmega))$	& Errore $L^2([0,T],H^1(\varOmega))$\\
		\midrule
		$4$				& $5.12828\cdot 10^{-1}$				& $1.91763$\\
		$5$				& $6.09235\cdot 10^{-2}$				& $8.40258\cdot 10^{-1}$\\
		$6$				& $2.00749\cdot 10^{-2}$				& $1.33768\cdot 10^{-1}$\\
		$7$				& $1.61626\cdot 10^{-3}$				& $8.43737\cdot 10^{-2}$\\
		\bottomrule
	\end{tabular}
	\caption{Norme degli errori al variare dei nodi nelle formule di quadratura per il calcolo dei polinomi.}
	\label{tab:grado}
\end{table}
Incidentalmente osserviamo che la scelta dei valori per il grado dei polinomi è stata piuttosto obbligata.
Infatti, per ragioni di tipo algoritmico, la libreria non può funzionare con meno di $4$ nodi di quadratura.
D'altra parte, una volta superati i $7$ nodi, le soluzioni sono talmente accurate da non permettere più un'analisi efficace, poiché intervengono errori di altro tipo a mascherare quelli che stiamo cercando.
In Figura~\ref{fig:poly} sono rappresentate le norme degli errori in corrispondenza dei diversi gradi di interpolazione.
\begin{figure}[tb]
	\centering
	\includegraphics[width=0.6\columnwidth]{poly}
	\caption{Logaritmi delle norme degli errori in norma $L^\infty([0,T],L^2(\varOmega))$ (in blu) e $L^2([0,T],H^1(\varOmega))$ (in rosso) al variare del numero di nodi nella formula di quadratura per il calcolo dei polinomi.}
	\label{fig:poly}
\end{figure}
In questo caso, per riuscire a mostrare con maggiore impatto la relazione tra le due variabili, abbiamo scelto di trasformare l'asse degli errori in scala logaritmica.
Secondo quanto predetto dalla teoria, la variazione dell'errore di interpolazione in funzione dell'ordine polinomiale è di tipo esponenziale.
In effetti, in figura abbiamo scelto di riportare le rette che meglio rappresentano l'andamento descritto dai punti e, come conseguenza, sembra proprio che i risultati teorici vengano validati, poiché c'è un buon accordo tra le rette e i punti.
Anche in questo caso ci siamo avvalsi del supporto di \texttt{R} e, dall'analisi diagnostica della regressione, risulta che le conclusioni appena descritte sono significative, nonostante la ridotta numerosità del campione.
D'altra parte, per noi non ha alcun interesse conoscere i veri coefficienti che descrivono le rette in modo analitico, poiché è sufficiente riconoscere l'andamento lineare in un piano logaritmico che, di conseguenza, mostra il \emph{trend} esponenziale citato poco prima.
Pertanto, possiamo concludere che aumentando il grado dei polinomi relativi alle formule di integrazione, entrambi gli errori sulla soluzione si riducono notevolmente, anche se l'errore in norma $L^\infty([0,T],L^2(\varOmega))$, come già succedeva nei casi precedentemente analizzati, si riduce di più.

Per concludere la nostra analisi, confrontiamo le prestazioni del $\vartheta$-metodo al variare del parametro di rilassamento.
Come dicevamo prima, poiché non esiste alcuna necessità di svolgere uno studio troppo sofisticato, focalizziamo la nostra attenzione sui $2$ metodi di maggiore interesse che siamo in grado di implementare per mezzo della nostra libreria:
\begin{itemize}
	\item il metodo di Crank-Nicolson, con $\vartheta=0.5$;
	\item il metodo di Eulero Implicito, con $\vartheta=1$.
\end{itemize}
Poiché il nostro obiettivo è quello di determinare l'ordine di convergenza degli schemi numerici, effettuiamo diverse simulazioni, variando il numero di passi necessari a risolvere nel tempo il problema numerico.
Per le successive analisi, in accordo con le condizioni utilizzate fin'ora, scegliamo di utilizzare una griglia spaziale composta di $16$ elementi spettrali per lato, con polinomi calcolati su $4$ nodi di integrazione.
Questa scelta ci permette, in pratica, di ridurre gli errori dovuti alla discretizzazione spaziale per mostrare la riduzione di quelli causati dallo schema di avanzamento in tempo.
I risultati di questa analisi sono sintetizzati in Tabella~\ref{tab:deltat}.
\begin{table}[tb]
	\centering
	\begin{tabular}{*4{c}}
		\toprule
		Metodo					& $\Delta t$	& Errore $L^\infty([0,T],L^2(\varOmega))$	& Errore $L^2([0,T],H^1(\varOmega))$\\
		\midrule
		\multirow{4}*{Crank-Nicolson}	& $0.1$	& $6.52108\cdot 10^{-2}$				& $9.33158\cdot 10^{-2}$\\
								& $0.05$	& $1.53864\cdot 10^{-2}$				& $2.31230\cdot 10^{-2}$\\
								& $0.025$	& $3.86838\cdot 10^{-3}$				& $7.61237\cdot 10^{-3}$\\
								& $0.0125$	& $9.635523\cdot 10^{-4}$				& $5.34705\cdot 10^{-3}$\\
		\midrule
		\multirow{4}*{Eulero Implicito}	& $0.1$	& $3.35351\cdot 10^{-1}$				& $5.50000\cdot 10^{-1}$\\
								& $0.05$	& $1.98667\cdot 10^{-1}$				& $3.00558\cdot 10^{-1}$\\
								& $0.025$	& $1.06809\cdot 10^{-1}$				& $1.57849\cdot 10^{-1}$\\
								& $0.0125$	& $5.54935\cdot 10^{-2}$				& $8.11585\cdot 10^{-2}$\\
		\bottomrule
	\end{tabular}
	\caption{Norme degli errori al variare del passo di avanzamento dello schema alle differenze finite.}
	\label{tab:deltat}
\end{table}
In ogni caso, per avere una rappresentazione più immediata dei risultati ottenuti, in Figura~\ref{fig:time} si possono osservare i logaritmi dei valori precedentemente inseriti in tabella.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Metodo di Crank-Nicolson}.]
		{\includegraphics[width=.45\columnwidth]{timeCN}} \quad
	\subfigure[\emph{Metodo di Eulero Implicito}.]
		{\includegraphics[width=.45\columnwidth]{timeEI}}
	\caption{Logaritmi degli errori di discretizzazione in funzione del logaritmo del passo di avanzamento dello schema alle differenze finite. In blu riportiamo i valori calcolati con la norma $L^\infty([0,T],L^2(\varOmega))$, mentre in rosso quelli relativi alla norma $L^2([0,T],H^1(\varOmega))$.}
	\label{fig:time}
\end{figure}
Dall'analisi delle immagini si può notare la presenza di una problematica già affrontata poco prima.
ll punto corrispondente all'errore di discretizzazione in norma $L^2([0,T],H^1(\varOmega))$ nel caso di Crank-Nicolson non sembra seguire l'andamento definito dalla teoria.
Come prima, la spiegazione di questo fenomeno risiede nella presenza degli altri errori, che attualmente non stiamo prendendo in considerazione.
In ogni caso, grazie agli strumenti regressivi già utilizzati, possiamo ricavare le rette che meglio approssimano i \emph{trend} definiti dai punti.
Nel caso dell'\emph{outlier} appena individuato, decidiamo di non considerarlo al momento della regressione, in modo da ottenere la retta che meglio si adatta ai dati definiti dagli altri $3$ punti.
In tutti gli altri casi le regressioni sono più che significative, nonostante la scarsa numerosità del campione, e ci permettono di trarre delle conclusioni importanti.
Infatti, il coefficiente angolare di queste rette ci consente di risalire all'ordine di convergenza degli schemi numerici alle differenze finite.
Una volta calcolati questi valori, li riportiamo in Tabella~\ref{tab:ordini}, in modo da verificare che le simulazioni numeriche siano in linea con i risultati teorici, i quali prevedono un andamento quadratico per il metodo di Crank-Nicolson, e uno lineare per tutti gli altri valori di $\vartheta$.
\begin{table}[tb]
	\centering
	\begin{tabular}{*3{c}}
		\toprule
		Metodo		& Ordine in norma $L^\infty([0,T],L^2(\varOmega))$	& Ordine in norma $L^2([0,T],H^1(\varOmega))$\\
		\midrule
		Crank-Nicolson	& $2.023380$							& $1.807852$\\
		Eulero Implicito	& $0.8681163$							& $0.921095$\\
		\bottomrule
	\end{tabular}
	\caption{Ordini di convergenza degli schemi numerici analizzati, rispetto al passo di avanzamento dello schema alle differenze finite.}
	\label{tab:ordini}
\end{table}
Dall'analisi dei valori in tabella, non possiamo che ritenerci soddisfatti degli ordini di convergenza, poiché ricalcano esattamente quelli che ci saremmo aspettati prima di intraprendere l'analisi.


\section{Un esempio applicativo}

Vediamo ora un'applicazione del nostro codice a un problema più complesso e concreto, ovviamente non risolvibile in maniera analitica.\\
Consideriamo il flusso di un gas in un tubo, all'interno del quale vengono posizionate delle resistenze che forniscono calore al gas stesso, traendo la propria energia da un impianto solare. 
Il tubo, in materiale metallico, è rivestito da una guaina isolante, che riduce la diffusione del calore nell'aria circostante.
Obiettivo dell'analisi è osservare la distribuzione del campo di temperatura nella sezione longitudinale del tubo.
Si suppone inoltre che, inizialmente, tutto l'apparato si trovi a temperatura ambiente, e che il flusso sia laminare, nonostante la velocità molto bassa.\\
Inquadrato in questo modo, si tratta di un problema parabolico di diffusione-trasporto, con condizioni al contorno di tipo Dirichlet, dato che consideriamo la temperatura ambiente imposta al bordo, a eccezione del lato di \emph{outflow}, dove si ha una condizione di Neumann omogenea (flusso nullo imposto).
Inoltre, intendiamo simulare l'effetto delle resistenze all'interno del tubo per mezzo di una forzante dipendente dal tempo.
\begin{equation}
	\left\{
	\begin{array}{ll}
		u_t - \mu \varDelta u + \vec{b}\cdot\nabla u = f & \qquad\textrm{in}\ \varOmega \times [0,T]\\[3mm]
		u(x,y,0) = u_0 & \qquad\textrm{in}\ \varOmega\\[3mm]
		u(x,y,t) = u_D & \qquad\textrm{su}\ \partial\varOmega\times [0,T]\\[3mm]
		\nabla u(x,y,t) \cdot \vec{n} = 0 & \qquad\textrm{su } \partial \varOmega\times [0,T].
	\end{array}
	\right.
\end{equation}
I parametri $\mu$ e $\vec{b}$, eventualmente dipendenti dallo spazio, rappresentano rispettivamente il coefficiente di diffusione e il coefficiente di trasporto, mentre $\vec{n}$ indica la normale uscente al dominio computazionale.

\subsection{Definizione del problema variazionale}

La definizione del problema è collocata nel file \texttt{tube.hpp}, contenuto nella cartella degli esempi.
Innanzitutto definiamo i parametri fisici legati al dominio.
\begin{c++}
double width(2);
double length(8);
double tube(1);
double thickness1(0.1);
double thickness2(0.2);
double thickness = thickness1 + thickness2;
double radius = tube - thickness;
double speed(3600);

const spectral::Point<double> LA(0., -width);
const spectral::Point<double> LB(length, width);
\end{c++}
Il dominio rettangolare ha una lunghezza orizzontale pari a \texttt{length} e una verticale pari a 2 volte \texttt{width}, dato che tutti i parametri verranno in seguito definiti tenendo presente la simmetria del problema. 
Lo spessore del tubo nel quale sono collocate le resistenze è pari a \texttt{thickness1}, mentre quello della guaina isolante è \texttt{thickness2}.
Infine, \texttt{speed} è un parametro di accelerazione, necessario per la riduzione dell'orizzonte temporale della simulazione.
In pratica, abbiamo bisogno di accelerare il processo, in modo che l'evoluzione del sistema sia totalmente osservabile simulando un breve intervallo temporale.
Di fatto, secondo questa trasformazione, $1$ secondo di simulazione corrisponde a \texttt{speed} secondi di processo fisico.
Nel caso specifico, il tempo che compare nei diversi parametri è misurato in ore anziché in secondi.\\
A questo punto è possibile definire la forzante tempo-dipendente:
\begin{c++}
struct TubeFunc {
    double operator()(double x,double y,double t) const {
        if (    (std::abs(y) > radius) && 
                (std::abs(y) < (tube-thickness2)) && 
                (((x > 0.2*length) && (x < 0.4*length)) || 
                ((x > 0.6*length) && (x < 0.8*length)))) {
            return 8*std::sin(2.*M_PI/24*t);
        } else {
            return 0;
        }
    }
};
\end{c++}
La forzante è nulla in tutto il dominio, ad eccezione di quattro zone collocate all'interno della parte metallica del tubo.
L'andamento è sinusoidale nel tempo, a rappresentare la dipendenza del riscaldamento dall'impianto solare, coerentemente con l'alternanza del giorno e della notte.\\ 
Dopodiché si prosegue con la definizione dei coefficienti di diffuzione e trasporto, anch'essi definiti come struttura, dato che non sono valori numerici costanti nel dominio.
\begin{c++}
struct Tubemu {
    double operator()(double x, double y) const {
        if (std::abs(y) > tube) {
            return 2.216e-5*speed;
        } else if (std::abs(y) > (radius+thickness1)){
            return 2.8e-7*speed;
        } else if (std::abs(y) > radius){
            return 1.12e-4*speed;
        } else { 
            return 2.204e-5*speed;
        }
    }
};
\end{c++}
Il coefficiente di diffusione assume un valore diverso a seconda del materiale a cui si riferisce.
Nel nostro caso abbiamo considerato i coefficienti relativi all'azoto per il gas all'interno del tubo, al rame per il tubo stesso, alla lana di vetro per la guaina isolante e all'aria per il gas che si trova nell'ambiente esterno al processo di diffusione del calore.
\begin{c++}
struct Tubeb1 {
    double operator()(double x, double y) const {
        if (std::abs(y) <= radius) {
            return -0.001*(y-radius)*(y+radius)*speed;
        } else {
            return 0;
        }
    }
};

struct Tubeb2 {
    double operator()(double x, double y) const {
        return 0;
    }
};
\end{c++}
Il coefficiente $\vec{b}$ è tale da simulare il trasporto dovuto a un flusso con profilo parabolico di velocità all'interno del fluido, unicamente nella direzione del moto del fluido stesso.
Per questo motivo, la componente lungo la direzione verticale è imposta pari a $0$.
Infine, come è naturale richiedere, all'esterno del tubo entrambe le componenti del campo di trasporto vengono scelte nulle.\\
A questo punto si può così definire la forma bilineare associata al problema, attraverso la classe \texttt{TubeBilinear}.
\begin{c++}
template<class MESH>
class TubeBilinear {
    public:
        TubeBilinear(MESH const & me): me(me) {
            for(int i=0; i<4; ++i) {
                bcond[i]=spectral::dirichlet;
            }
                bcond[1]=spectral::neumann;
        }
        // Boundary conditions
        spectral::bctype bcond[4];

        struct neumann{
            double operator()(double const & x, 
                double const & y) const {
                return 0;
            }
        };
        struct dirichlet{
            double operator()(double const & x, 
                double const & y) const {
                return 300;
            }
        };
\end{c++}
Le condizioni al contorno vengono assegnate attraverso le strutture \texttt{on}, sia per il tipo Neumann sia per il tipo Dirichlet.
Grazie alla templetizzazione del codice, è possibile definire il metodo per la condizione di Dirichlet per ogni lato, mentre per quella di Neumann è sufficiente assegnare il numero del lato alla variabile intera.
\begin{c++}
template<unsigned int N, typename T=void>
struct on : public dirichlet {};

template<typename T>
struct on<1,T> : public neumann {};
\end{c++}
In realtà, la forma bilineare vera e propria viene costruita a partire dall'operatore~\texttt{()}, il quale contiene le informazioni legate ai termini di diffusione e trasporto, secondo la formulazione debole.
Viene imposta anche la condizione iniziale, che assegna la temperatura di $300$~K a tutto il dominio.
Infine, per completezza riportiamo anche gli elementi privati della classe \texttt{TubeBilinear}.
\begin{c++}
template<class U,class V> inline
double operator()(U const & u,V const & v,int n) const {
    return
        + me.prodL2(me.restriction(mu,n),me.Dx(u,n),me.Dx(v,n),n)
        + me.prodL2(me.restriction(mu,n),me.Dy(u,n),me.Dy(v,n),n)
        + me.prodL2(me.restriction(b1,n),me.Dx(u,n),v,n)
        + me.prodL2(me.restriction(b2,n),me.Dy(u,n),v,n);
}
struct initial{
    double operator()(double const & x, double const & y) const {
        return 300;
    }
};
private:
    MESH const & me;
    Tubemu mu;
    Tubeb1 b1;
    Tubeb2 b2;
};
\end{c++}
La macro \texttt{TUBE\_PROBLEM} contiene tutte le informazioni relative al problema. 
Dato che in questo caso si utilizza una forzante tempo dipendente, è sufficiente assegnare il valore \texttt{true} all'ultimo flag della classe \texttt{Functional}, anziché definirne una nuova.
\begin{c++}
#define TUBE_PROBLEM TubeBilinear<spectral::Mesh>,
		Functional<spectral::Mesh,TubeFunc,true>
\end{c++}
Arrivati a questo punto, possiamo affermare di aver creato tutte le strutture necessarie per la definizione del problema e, finalmente, possiamo concentrarci sulla sua risoluzione numerica.

\subsection{Definizione dei parametri e risoluzione numerica}

Possiamo procedere con la costruzione del \texttt{main()}, incluso nel file \texttt{tube.cpp}, che conterrà tutte le istruzioni relative all'impostazione dello schema numerico risolutivo.
\begin{c++}
// Number of elements along the x-axis
int nx(50);
// Number of elements along the y-axis
int ny(50);
// Number of nodes
int nnodes(4);
// Time step
double deltat(0.5);
// Number of time steps
int ntime(96);
// Relaxation parameter
double theta(1);
// nmax for the newton method
int newton_nmax((int)1e4);
// Order of the problem in the week form
int order(1);
// Tolerance parameters for the newton method
double legDelta(1e-15);
double legEps(1e-15);
\end{c++}
Abbiamo scelto di assegnare $50$ elementi su ogni lato, e di utilizzare un passo di discretizzazione temporale pari a $0.5$~ore, per un totale di $96$ iterazioni, in modo da ottenere un orizzonte temporale di $48$~ore.
In questo modo la forzante risulterà attiva nelle ore diurne, con un andamento sinusoidale che intende simulare l'effetto della variazione dell'intensità dell'energia solare nel corso della giornata.\\
Una volta che tutti i parametri sono stati definiti, è possibile passare alla risoluzione del problema vera e propria.
Per prima cosa, il problema viene inizializzato.
\begin{c++}
Problem<TUBE_PROBLEM> pb(LA,LB,nx,ny,nnodes,deltat,
                             theta,newton_nmax,
                             legDelta,legEps,order);
		pb.initialize();
\end{c++}
Diversamente dal caso test nell'esempio precedente, in cui avevamo la soluzione esatta del problema, in questo caso non ci è possibile calcolare le norme degli errori.
In realtà non è nemmeno l'obiettivo di questa simulazione che, di fatto, ha un interesse puramente applicativo e dimostrativo, nei confronti delle potenzialità della libreria.\\
Procediamo, pertanto, con la risoluzione del problema evolutivo.
\begin{c++}
double actualtime(0);
for (int t=0; t<ntime; ++t) {
    actualtime = deltat*(t+1);
    pb.solve(actualtime);
}
\end{c++}
Dato che la forzante è dipendente dal tempo, il suo contributo deve essere valutato ad ogni istante temporale, anche se le matrici del problema non cambiano la loro struttura.
Per tenere conto di questo fatto, è necessario passare in ingresso al metodo \texttt{solve} l'istante temporale in cui si va a calcolare la soluzione.
Questo stratagemma, grazie allo stile di programmazione che è stato adottato in fase di scrittura della libreria, ci consente di utilizzare un metodo con lo stesso nome di quello che avremmo usato nel caso di forzante stazionaria nel tempo.
In effetti, i due metodi hanno lo stesso nome ma diversa firma, poiché quello che usiamo ora prende in ingresso una variabile di tipo \texttt{double}.
Questo vuole essere un punto di forza del codice e, indubbiamente, mette in luce la versatilità che la libreria può dimostrare nelle situazioni più disparate.

\subsection{Risultati delle simulazioni}

I risultati della risoluzione numerica possono essere rappresentati con semplicità, sfruttando gli script automatici generati dal codice. 
Il video che si ottiene a partire dallo script \texttt{heat2mencoder.sh} visualizza l'intera evoluzione nel tempo del campo di temperatura, all'interno del dominio preso in considerazione.
Lo script \texttt{heat2gnuplot.sh}, invece, permette la visualizzazione della soluzione a un certo istante temporale, grazie alla sintassi da terminale
\begin{verbatim}
chmod 744 heat2gnuplot.sh
./heat2gnuplot.sh N
\end{verbatim}
dove \texttt{N} è il numero dell'iterazione che l'utente intende rappresentare in un'immagine, la quale sarà generata nel formato vettoriale \texttt{pdf}.
\begin{figure}[tb]
	\centering
	\subfigure[\emph{Istante con forzante massima}.]
		{\includegraphics[width=\columnwidth]{tubeSolMax}} 
	\subfigure[\emph{Istante con forzante minima}.]
		{\includegraphics[width=\columnwidth]{tubeSolMin}}
	\caption{Immagini della soluzione in presenza della guaina isolante, a diversi istanti temporali.}
	\label{fig:tubeSol}
\end{figure}
Grazie alle istruzioni appena mostrate, siamo in grado di rappresentare in Figura~\ref{fig:tubeSol} il campo di temperatura, in tutto il dominio computazionale, relativo all'istante in cui la forzante è massima e a quello in cui è minima.
Si osserva come il calore, che è generato dalle resistenze nel tubo metallico, viene trasferito per diffusione al gas all'interno del tubo, mentre la guaina isolante evita l'innalzamento della temperatura dell'aria all'esterno.
Inoltre, in entrambe le situazioni, è evidente l'effetto del trasporto all'interno del tubo: anche quando le resistenze non sono più attive, l'azoto trasportato all'uscita della condotta risulta decisamente surriscaldato.
\begin{figure}[tb]
        \centering
        \subfigure[\emph{Istante con forzante massima}.]
                {\includegraphics[width=\columnwidth]{tubeSolMaxNo}} 
        \subfigure[\emph{Istante con forzante minima}.]
                {\includegraphics[width=\columnwidth]{tubeSolMinNo}}
        \caption{Immagini della soluzione in assenza della guaina isolante, a diversi istanti temporali.}
        \label{fig:tubeSolNo}
\end{figure}
Invece, nel caso in cui la guaina non è presente, come si può osservare in Figura~\ref{fig:tubeSolNo}, il comportamento globale è differente.
La diffusione coinvolge anche l'aria all'esterno del tubo, che però non è soggetta al processo di trasporto, con un conseguente surriscaldamento della stessa.
Inoltre, i valori di temperatura del gas interno risultano diminuire rispetto al caso precedente.
In conclusione, possiamo affermare che la presenza di una guaina isolante è necessaria se l'intento del progettista è quello di simulare un impianto di riscaldamento ad energia solare.\\
Per finire, dobbiamo necessariamente sottolineare che si tratta comunque di un'approssimazione grossolana di un impianto reale.
Ad esempio, la velocità del gas all'interno della condotta è stata ridotta in maniera artificiosa, in modo da riuscire a mostrare in contemporanea gli effetti della diffusione e del trasporto, difficilmente presenti nello stesso momento se si utilizzano parametri realistici.
In ogni caso, lo scopo principale di questa simulazione è soltanto quello di fornire un esempio per la realizzazione di altri problemi tempo-dipendenti, sfruttando la libreria \texttt{Heat-Spectral}, nonché quello di mostrare la flessibilità del codice, in grado di trattare, ad esempio, condizioni al contorno miste, omogenee o non omogenee, forzanti costanti o variabili nello spazio e nel tempo, oltre a coefficienti costanti o variabili in spazio.
